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ABSTRACT 

We present improved synthesis models of the evolving spectrum of the UV/X-ray diffuse back- 
ground, updating and extending our previous results. Five new main components are added to our 
radiative transfer code CUBA: (1) the sawtooth modulation of the background intensity from res- 
onant line absorption in the Lyman scries of cosmic hydrogen and helium; (2) the X-ray emission 
from the obscured and unobscurcd quasars that gives origin to the X-ray background; (3) a piccewise 
parameterization of the distribution in redshift and column density of intergalactic absorbers that fits 
recent measurements of the mean free path of 1 ryd photons; (4) an accurate treatment of the pho- 
toionization structure of absorbers, which enters in the calculation of the helium continuum opacity 
and recombination emissivity; and (5) the UV emission from star-forming galaxies at all redshifts. 
We provide tables of the predicted H i and He ii photoionization and photolieating rates for use, 
e.g., in cosmological hydrodynamics simulations of the Lya forest, and a new metallicity-dependent 
calibration to the UV luminosity density-star formation rate density relation. A "minimal cosmic 
reionization model" is also presented in which the galaxy UV emissivity traces recent determinations 
of the cosmic history of star formation, the luminosity-weighted escape fraction of hydrogen-ionizing 
radiation increases rapidly with lookback time, the clumping factor of the high-redshift intergalactic 
medium evolves following the results of hydrodynamic simulations, and Population HI stars and mini- 
quasars make a negligible contribution to the metagalactic flux. The model provides a good fit to the 
hydrogen-ionization rates inferred from fiux decrement and proximity effect measurements, predicts 
that cosmological H ii (He iii) regions overlap at redshift 6.7 (2.8), and yields an optical depth to 
Thomson scattering, Tos = 0.084 that is in agreement with WMAP results. Our new background 
intensities and spectra are sensitive to a number of poorly determined input parameters and suffer 
from various degeneracies. Their predictive power should be constantly tested against new observa- 
tions. We are therefore m aking our redshift-dependcnt UV/X emissi vities and CUBA outputs freely 



available for public use at littp://www. ucolick.org/~pmadau/CUBA 
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1. INTRODUCTION 



The reionization of the all-pervading intergalactic medium (IGM) is a landmark event in the history of cosmological 
structure formation. Studies of Gunn-Peterson absorption in the spectra of distant quasars show that hydrogen is 
highly photoionized out to redshift z > 6 (e.g., Fan, Carilli, & Keating 2006a; Songaila 2004), while polarization data 
from the Wilkinson Microwave Anisotropy Probe (WMAP) constrain the redshift of a sudden reionization event to be 
significantly higher, z = 10.5 ± 1.2 (Jarosik et al. 2011). It is generally thought that the IGM is kept ionized by the 
integrated UV emission from active nuclei and star-forming galaxies, but the relative contributions of these sources as 
a function of epoch are poorly known. Because of the high ionization threshold (54.4 eV) and small photoionization 
cross section of He ii, and of the rapid recombination rate of He iii, the double ionization of helium is expected to be 
completed by hard UV-cmitting quasars around the peak of their activity at z w 2.5 (e.g., Madau & Meiksin 1994; 
Sokasian, Abel, & Hernquist 2002; McQuinn et al. 2009), much later than the reionization of H i and He i. At z > 3, 
the declining population of bright quasars appears to make an increasingly small contribution to the 1 ryd radiation 
background, and it is believed that massive stars in galactic and subgalactic systems may provide the additional 
ionizing fiux needed at early times (e.g., Madau, Haardt, & Rees 1999; Gnedin 2000; Haehnelt et al. 2001; Wyitlie & 
Loeb 2003; Meiksin 2005; Trac & Cen 2007; Faucher-Giguere et al. 2008a; Gilmore et al. 2009; Robertson et al. 2010). 
This idea may be supported by the detection of escaping ionizing radiation from individual Lyman-break galaxies at 
z - 3 (e.g., Shapley et al. 2006). 

Despite much recent progress, a coherent description of the thermal state and ionization degree of the IGM remains 
elusive. The intensity and spectrum of the cosmic ultraviolet background remain one of the most uncertain yet 
critically important astrophysical input parameters for cosmological simulations of the IGM and for interpreting 
quasar absorption-line data and derive information on the distribution of primordial baryons (traced by H i. He i. 
He II transitions) and of the nucleosynthetic products of star formation (C iii, C iv, Si iii, Si iv, O vi, etc.). This is 
the fourth paper in a series aimed at a detailed study of the generation and reprocessing of photoionizing radiation 
in a clumpy universe, and of the transfer of energy from this diffuse background flux to the IGM. In Paper I (Madau 
1995) we showed how the stochastic attenuation produced by neutral hydrogen along the line of sight affects the colors 
of distant galaxies. In Paper II (Haardt & Madau 1996) we developed CUBA, a radiative transfer code that followed 
the propagation of Lyman-continuum (LyC) photons through a partially ionized inhomogeneous IGM. CUBA outputs 
have been extensively used to model the Lya forest in large cosmological simulations (e.g., Tytler et al. 2004; Theuns 
et al. 1998; Dave et al. 1997; Zhang et al. 1997). In Paper III (Madau et al. 1999) we focused on the candidate 
sources of photoionization at early times and on the history of the transition from a neutral IGM to one that is 
almost fully ionized. In this paper we describe a new version of CUBA and use it to compute improved synthesis 
models of the UV/X-ray cosmic background spectrum and evolution, combining, updating, and extending many of our 
previous results in this field. The five main upgrade to CUBA are: (1) the sawtooth modulation from resonant line 
absorption in the Lyman series of intergalactic helium as well as hydrogen; (2) the X-ray emissivity from the obscured 
and unobscured populations of active galactic nuclei (AGNs) that gives origin to the X-ray background; (3) an up-to- 
date piecewise parameterization of the distribution in column density of intervening absorbers, which establishes the 
"super Lyman-limit systems" as the dominant contributors to the hydrogen LyC intergalactic opacity; (4) an accurate 
treatment of the absorber photoionization structure, entering in the calculation of the helium continuum opacity and 
recombination emissivity of the clumpy IGM; and (5) the UV flux from star-forming galaxies at all redshifts. 

The plan is as follows. In § 2 we review the basic theory of cosmological radiative transfer in a clumpy universe. 
§ 3 and § 4 discuss the distribution of absorbers along the line of sight and their photoionization structure. The 
recombination radiation from the clumpy IGM is calculated in §5. In §6 and §7 we compute the UV and X-ray 
emissivity from quasars, and in § 8 the UV emissivity from star-forming galaxies. An overview of the main results 
generated by the updated CUBA radiative transfer code is given in § 9. Finally, we summarize our findings in § 10. 
Unless otherwise stated, all results shown below will assume a (f^M, ^^A, ^b, h) = (0.3, 0.7, 0.045, 0.7) cosmology. Note 
that, while the source volume emissivities must be evaluated in a given cosmological model, the resulting background 
intensity does not explicitly depend on the choice of cosmological parameters. 



We start by summarizing the basic theory describing the propagation of ionizing radiation in a clumpy, primordial 
IGM (e.g.. Paper I; Paper II; Madau & Haardt 2009). The equation of cosmological radiative transfer describing the 
time evolution of the space- and angle-averaged monochromatic intensity is 



where H{z) is the Hubble parameter, c the speed of the light, is the absorption coefficient, and the proper volume 
emissivity. The integration of equation (|^) gives the background intensity at the observed frequency i^o, as seen by an 
observer at redshift Zg, 



where v = Vo{l + z)/{l + Zq), \dt/dz\ = H ^(1 4- z) ^, r = — ln(e '^) is the effective absorption optical depth of a 
clumpy IGM, and is the proper volume emissivity. 



2. COSMOLOGICAL RADIATIVE TRANSFER 




(1) 




(2) 
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2.1. Photoelectric absorption 
In the case of LyC absorption by Poisson-distributed systems, the effective opacity between Zq and z is 

f,K,Zo,z)= / dz' dNmf{Nm,z'){l~e-^^), (3) 

JO 

where /(iVni, z') is the bivariate distribution of absorbers in redshift and column density along the line of sight, Tc is 
the continuum optical depth at frequency v' = 1^0(1 + z')/{l + z^) through an individual absorber, 

T^{i^') = Niiia^iiiu') + A^HoiCTHci(i^') + N^cuCTHcuii^'), (4) 
where and are the column densities and photoionization cross sections of ion i. 

2.2. Resonant absorption 

Besides photoelectric absorption, resonant absorption by the hydrogen and helium Lyman series will produce a 
sawtooth modulation of the radiation spectrum (Madau & Haardt 2009; Haiman, Rees, & Loeb 1997). Continuum 
photons that are redshifted through the Lya frequency, Va, are resonantly scattered until they redshift out of resonance: 
the only two Lya line destruction mechanisms, two-photon decay and O iii Bowen fluorescence (Kallman & McCray 
1980), can typically be neglected in the low metallicity, low density IGM. This is not true, however, for photons passing 
through a higher order Lyman-series resonance, which will be absorbed and degraded via a radiative cascade rather 
than escaping by redshifting across the line width. Since the line absorption cross section is a narrow, strongly peaked 
function, the effective line absorption optical depth for a photon observed at (zoj^^o < ^^n) that passed through a 
resonance at redshift z„ = (1 + Zo){vn/^'o) — 1, can be written as 

T„(z„) - (1 + z„)— / dNmf{Nm,Zn)Wn, (5) 
c Jo 

where = i/q, x 4(1 — n~^)/3 is the frequency of the Is — t- np Lyman-series transition {n > 3) and Wn is the rest 
equivalent width of the line expressed in wavelength units. This opacity is dominated by systems having line center 
optical depths of order unity, i.e., which lie at the transition between the linear and the flat part of the curve of growth. 

Consider, for example, radiation observed at frequency below the Ly/3 of hydrogen or helium, Vo < vp. Photons 
emitted between Zo and zp ~ [1 + Zo)(i'i3/i'o) — 1 can reach the observer without undergoing resonant absorption. 
Photons emitted between zp and = (1 + Zo){i'~f / Vo) ~ 1 pass instead through the Ly/3 resonance at and are 
absorbed. Photons emitted between z^ and 25 = (1 + Zo){i's/i^o) — 1 pass through both the Ly/3 and the Ly7 
resonances before reaching the observer. The background intensity can then be written as (Madau & Haardt 2009) 

00 

J^^zo) ^ 0{zo,zp) +0{zfi,z^)e~-^^ +0{z^,zs)e-^^~^-- +.... + 0(zz„ 00) exp(- ^ f„), (6) 

n— 3 



where we have denoted with the symbol 0{zi,Zj) the "Gibers' integrals" on the right hand side of equation 




calculated between redshifts Zi and Zj > Zi and with equal to the relevant continuum opacity, 

0(z„ z,) ^ ^ \dt/dz\dz^-j^±^e.{z)e~'^. (7) 

In equation (^, = (1 + Zo){vlIvo) — 1, vl is the frequency at the Lyman limit, and fp,f-y,fs, ... are the Lyman- 
series effective opacities at redshifts zp, z-y, zs, .... In the case of resonant absorption by H i, the LyC optical depth 
fc is zero in all O-integrals except the last, while in the case of He 11 all terms must include photoelectric absorption 
by H I and He i (as well by He 11 in the last term). Equation (||) is easily generalized to higher observed frequencies, 
Vn < < to read 

00 k 

JuAzo) ^ 0{zo,Zn+i) + ^ 0(zfc,Zfc+i) exp(- ^ f;). (8) 

k—n+1 l—n+l 

Note how, in the case of large resonant opacities, only sources between the observer and the "screen" redshift z„ = 
(1 -|- Zo){vn/vo) — 1 corresponding to the frequency of the nearest Lyman-series line above Vo will not be blocked from 
view: the background energy spectrum will show a series of discontinuities, peaking at frequencies just above each 
resonance, as the first integral in equation (|^) extends over the largest redshift path, and going to zero at resonance. 

2.3. Lyman series cascades 

Each photon absorbed through a Lyman series resonance causes a radiative cascade that ultimately terminates either 
in a Lya photon or in two-photon 2s — >• Is continuum decay. In the former case the photon scatters until it is redshifted 
out of resonance, in the latter the photons escape to infinity without further interactions. Consider, for example, the 
absorption of a Ly/3 photon. The excited 3p level is depopulated via 3p 2s decay (Ha). In the low-density IGM, 
coUisional Z-mixing of the 2s — 2p levels (Seaton 1959) is negligible, and the cascade can only terminate in two-photon 
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2s — > Is emission (Hirata 2006). Without /-mixing, the quantum selection rules forbid a Ly/3 photons from being 
converted into Lya: by contrast, excitation of the 4p level by absorption of Ly-f can decay via the 3s or 3d levels to 
the 2p and ultimately produce Lya. More generally, the fraction, /„, of decays from an np state that generates Lya 
photons can be determined from the selection rules and the decay probabilities. This fraction is found to increase as 
/„ = (0, 0.2609, 0.3078, 0.3259, ...) for n = (3, 4, 5, 6, ...), and to asymptote to 0.359 at large n (Pritchard & Furlanetto 
2006). Note that this is valid in the approximation that the IGM is optically thick to higher-order Lyman-series 
transitions. 

What is the Lya diffuse flux produced by these Lyman-series cascades? Let Ju„[z) be the background intensity 
measured just above the H i or He ii Lyn resonance at redshift z. The Lyn flux that is absorbed and converted into 
Lya is fnJy^{z)[l — e""^"*^^-*], and the proper Lya volume cmissivity generated by this process can be written as 

e-[z) = J.„(z)[l - e-^-f^)] ^^^^^ \dz/dt\, (9) 

where 5{x) is the Dirac delta function. The additional flux observed at frequency Vq < Va and redshift Zq from this 
process is then 

^J:S^o)={j-^ e-^^^''°'^°'^"'{/«J.„(^a)[l-e^^"(^")]} (10) 

(Madau & Haardt 2009), where 1 -I- Zq, = (1 4- Zo){va/vo) and the LyC optical depth is zero in the case of resonant 
absorption by H I. When summing up over all Lyman series lines, the term in square brackets must be replaced by 
^j^^g{/„Ji,^(zQ,)[l — e~'^"(^°)]}. The same Lyn cascade also produces a two-photon continuum with emissivity given 

by 

e:(z) = -(1 - /„)J.„(z)[l - e-^-f^)] \dz/dtl (11) 

C (1 + 2:) 

where the two-photon emission function is expressed in photons per unit frequency interval and is symmetric about 



We note that the underlying assumption in equations (|9|) and (11) is that every absorber is a source of unprocessed 
Lya line and two-photon continuum radiation, i.e., that these photons escape into intergalactic space without ap- 
preciable local absorption. In the case of He ii Lya emission, this requires negligible "in situ" destruction from dust, 
metals (O iii Bowen fluorescence), and photoelectric absorption by H i, so that the Lya photons diffuses into the wings 
and eventually escape from the production site into the IGM. This is a good approximation at the low metallicities 
that characterize intergalactic absorbers (Kallman & McCray 1980), even more so since the reprocessing of Lyman 
series photons occurs in a "skin" layer at the surface of an absorption system. In § 4 we will show that this is a poor 
approximation in the case of the reprocessing of LyC radiation, a proper treatment of which requires a numerical 
solution of the radiative transfer equation within individual absorbers. 



TABLE 1 

Parameters of the distribution of intergalactic absorbers 



Absorbers class 


log (JVHi/cm-2) 


P 


A [cm-2(/3-i)] 


7 


redshift 


Lya forest 


11 - 15 


1.5 


107.079 


3.0 


1.56 < < 5.5 




11 - 15 


1.5 


]^q8.238 


0.16 


z < 1.56 




11 - 15 


1.5 


101.470 


9.9 


2 > 5.5 




15 - 17.5 


2.0 


1014-580 


3.0 


1.56 < 2 < 5.5 




15 - 17.5 


2.0 


1015.740 


0.16 


2 < 1.56 




15 - 17.5 


2.0 


108.970 


9.9 


2 > 5.5 


LLSs 


17.5 - 19 










SLLSs 


19 - 20.3 


1.05 


10-0.347 


1.27 


2 > 1.56 




19 - 20.3 


1.05 


100.107 


0.16 


2 < 1.56 


DLAs 


20.3 - 21.55 


2.0 


1018.940 


1.27 


2 > 1.56 




20.3 - 21.55 


2.0 


1019.393 


0.16 


2 < 1.56 



3. DISTRIBUTION OF ABSORBERS ALONG THE LINE OF SIGHT 

The effective opacity of the IGM has traditionally been one of the main uncertainties affecting calculations of the 
UV background. Our improved model uses a piecewise power-law parameterization for the distribution of absorbers 
along the line of sight, 

f{Nui,z)^AN^fil + z)\ (12) 
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Fig. 1. — Left: The assumed distribution of absorbers in H I column density at redshift z = 2 {bottom curve), z = 3.5 (middle curve), 
and 2 = 5 {top curve). For clarity, we have multiplied the top and bottom curves by 50 and 1/50, respectively. Right: The quantity 
-'^hi/C-'Vhi, z)[l — exp(— AfHio"9i2] at the same redshifts, showring the dominant contribution of the optically thick LLSs and SLLSs to the 
intergalactic opacity at 1 ryd. 



and is designed to reproduce accurately a number of recent observations: 

• Over the column density range 10^^ < Nm < 10^^ cm~^, we use (A,/3,7) = (1.2 X 10^1.5,3.0), where the 
normalization A is expressed in units of cm^^('^^^\ and (3 = 1.5 is chosen following, e.g., Tytler (1987). As noted, 
e.g., by Meiksin & Madau (1993), Petitjean et al. (1993), and Kim et al. (1997, 2002), /(A^m) starts to steepen 
from the empirical —1.5 power law at A'hi > lO^^'^cm"^. Here, we assume (A,/?, 7) = (3.8 x 10^"', 2.0, 3.0) 
for 10^^ < A'hi < 10^^'^ cm~^. A "curve of growth" analysis (providing the relationship between equivalent 
width and column density) with Doppler parameter & = 32 kms~^, together with equation (^ and the above 
distribution of Lya-forest clouds, produces a Lya effective opacity fa = 0.0015(1 + z)^, in agreement with the 
best-fits of Fauchcr-Giguere et al. (2008b) after metal correction. 

• At the other end of the column density distribution, a recent survey of damped Lya systems (DLAs) by Prochaska 
& Wolfe (2009) (see also Guimaraes et al. 2009) yields dN/dz = J dNmfiNm,z) = 0.294 DLAs per unit 
redshift at (z) = 3.7 above A^hi = lO^*^ "^ cm~^. With a power-law exponent /3 = 2 down to a break column of 
A^Hi = lo^i-ss cm~^ (Prochaska & Wolfe 2009), and with an incidence per unit redshift oc (1 + z)^-^^ (Rao et al. 
2006), the parameters for the DLAs becomes (A,/3,7) = (8.7 x 10^^,2,1.27). 

• For absorbers with 10^^ < A^hi < 10^°'^ cm~^ (the so-called "super Lyman-limit systems", or SLLSs), we use 
O'Meara et al. (2007), who find dN/dz = 0.97 SLLSs per unit redshift at (z) = 3.5 above Nm = lO^^ cm'^. 
Matching with the DLAs abundance then requires {A, (3, 7) = (0.45, 1.05, 1.27) for the SLLSs. 

• There is obviously a significant mismatch between the power-law exponent for the Lya clouds (7 = 3) and the 
SLLSs (7 = 1.27). Continuity then requires the shape of f{Nm,z) to change with redshift over the column 
density range of the Lyman-limit systems (LLSs), 10^^'^ < A^hi < 10"'^^ cm~^. In this interval of column densities 
we match the distribution function with a power law of redshift-dependent slope. The procedure yields the slopes 
/? = 0.47,0.61,0.72,0.82 at redshifts z = 2,3,4,5, respectively, in agreement with Prochaska et al. (2010) who 
find for the LSSs (3 = Q.Sto.t at z « 3.5. 

• The ensuing f{N-iii,z) distribution is shown in the left panel of Figure for z — 2,3.5,5 where, for clarity, we 
have multiplied the values at the highest and lowest redshift by 50 and 1/50, respectively. Its shape is similar 
to the distribution inferred by Prochaska et al. (2010). In the right panel of the same figure we have plotted the 
quantity Niiidfc/{dzdNiii)\u=ugi2 = -^hi/(A'hi, ^)[1 — cxp(— A^hiO'912], i-e., the effective optical depth at 1 ryd 
per unit redshift per unit logarithmic interval of hydrogen column. This shows the dominant contribution of the 
LLSs and SLLSs to the LyC opacity. 
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• The above paramcterizations reproduce well the observations at 2 < z < 5. At low redshift, however, Hubble 
Space Telescope (HST) data show that the forest undergoes a much slower evolution. Following Weymann et al. 
(1998) we take 7 = 0.16 in the interval < z < ziow and dN/dz = 34.7 at z = above an equivalent width of 0.24 
A (corresponding to a column of lO^'^-^'^ cm~^ for b = 32 kms~^). We derive (A, /3,7) = (1.73 x 10^,1.5,0.16) 
for 10" < A^Hi < 10^^ cm-2 and (A,/3,7) = (5.49 x 10l^2,0.16) for lO^^ < iVm < lO"-^ cm^^ at all redshifts 
below ziow = 1.56. We use a broken power-law for the redshift distribution of the SLLSs and DLAs as well; 
assuming that the same 7 = 0.16 slope and transition redshift ziow inferred for the forest also hold in the case of 
the thicker absorbers, we derive a normalization at z < z\o^ oi A = 1.28 for the SLLSs and A = 2.47 x 10^^ for 
the DLAs. This yields dN/dz = 0.74 absorbers above Nm = lO^^'^ cm'^ at (z) = 0.69, in agreement with the 
value measured by Stengler-Larrea et al. (1995), dN/dz ~ 0.70 ± 0.2. 

• Above z = 5.5, the spectra of the highest redshift quasars known show an accelerated evolution in the Lya 
opacity of the IGM, Tq = 2.68[(1 + z)/6.5]^°'^ (Fan et al. 2006b), indicating a sharp increase in in the average 
neutrality of the universe. This can be mimicked by assuming for the forest the values [A, /3, 7) = (29.5, 1.5, 9.9) 
(10" < Nm < 10^^ cm-2) and (A,/3,7) = (9.35 x 10^2,9.9) (lO^^ < Nm < lO^"^'^ cm'^) above redshift 5.5. 




redshift redshift 



Fig. 2. — Left: The predicted proper mean free path at 1 ryd {solid line) together with the measurements of Prochaska et al. (2009) 
(crosses). Right: Evolution of the observed effective Lya optical depth, — ln(T), where T is the transmitted flux ratio. Data points are 
from Schaye et al. (2003; red filled circles), Songaila (2004; magenta empty squares), Faucher-Giguere et al. J^OOSb; green filled tefuares), 
and Fan et al. (2006b; blue filled circles). The solid line shows the Lya opacity, fa, predicted by equations (jsh (for n = 2) and (112), and 
using a curve-of-growth analysis corresponding to a Dopplcr parameter b = 32 kms~^. 



The parameters of the adopted distribution of intergalactic absorbers are summarized in Table 1. 

3.1. Mean free path of hydrogen-ionizing radiation 
Inserting our f{Nm, z) in equation (^), we can compute the (proper) LyC mean free path for 1 ryd photons as 

dz 

X'^l = c\dt/dz\x—\,=,,,,. (13) 

This is plotted in the left panel of Figure || in the redshift range 3.5-4.5. At z = 3.5, the major contributors to the 
LyC opacity are, in order of decreasing magnitude, the high column-density Lya forest (10^^ < A'hi < 10^^'^ cm~^, 
32%), the SLLSs (28%), the LLSs (20%), the low column-density Lya forest (iVm < 10^^ cm-^, 12%), and the DLAs 
(8%). A new method to directly measure the IGM LyC opacity along quasar sight lines has been recently presented 
by Prochaska, Worseck, & O'Meara (2009). The approach analyzes the "stacked" spectrum of 1,800 quasars drawn 
from the Sloan Digital Sky Survey (SDSS) to give an empirical determination of the mean free path A^^^jp. Our new 
opacity model agrees very well with the measurements of Prochaska et al. (2009), and produces a continuum opacity 
that is approximately half of that adopted in Paper II. The right panel of the same figure shows how our model also 
provides a good fit to the Lya quasar transmission data over the entire redshift range 2 < z < 6. 

For a single population of absorbers described by equation (p^ ) , the mean free path scales with frequency and redshift 
as Amfp(i^, z) oc (i//i/gi2)^^'^~"'^'iJ~^/(l -l- z)'^^^ . Givcn the multi-component distribution summarized in Table 1, we 
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can readily compute the mean free path of ionizing radiation in the range 13.6 < hv < 48.4 cV under the assumption 
that He i continuum absorption at energies above 24.2 cV can be neglected (photons between 48.4 and 54.4 eV are 



reprocessed by He ii Lyman series resonance absorption, see § 2.2). For ease of use in analytical calculations, we fit 
our numerical results for the mean free path as 

Amtp(i^, z) = c\dt/dz\/^z = c\dt/dz\A{s){l + z)'^^"' (14) 

where s = v/vqi2- Both the normalization A{s) and the exponent 7(5) are well fit by third order polynomials of the 
form 

[A{s),i{s)] =P3(s'-l)+P2(s'-l)+Pi(s-l)+Po. (15) 

Numerical values of the best-fit polynomial coefficients are given in Table 2: the fitting function is adjusted to be 
continuous in value at the redshifts where it changes slope. As discussed above, the fit is only valid in the range 
1 < s < 3.56. Close to the hydrogen Lyman edge, and at early enough epochs, only "local" radiation sources - 
sources within a mean free path of a few tens of Mpc - contribute to the ionizing background intensity, and one can 
neglect cosmological effects such as source evolution and frequency shifts. In this "source-function" approximation, 
47rJ9i2(z) « €912 (z) A^,if2 (z). 



TABLE 2 

Fitting parameters for the hydrogen LyC mean free path 



parameter pa p2 Pi Po 

A 0.0509 -0.406 1.167 1.076 
< 2 < 1.56 7 0. 0. 0. -0.160 

1.56 < 2 < 5.5 7 0.0593 -0.519 1.586 -2.104 
z > 5.5 7 0.122 -1.356 5.998 -8.423 



4. PHOTOIONIZATION STRUCTURE OF ABSORPTION SYSTEMS 

The ionization state of individual absorbers enters in calculations of the He i and He 11 opacities and of the continuum 
and line recombination radiation from hydrogen and helium. Under the assumption of photoionization equilibrium 
(generally accurate for quasar absorbers, see Paper II), in a pure H/He gas illuminated by a local radiation intensity 
J7y, the ion fractions Ihi, ^HcI, and IhcII can be written in implicit form as 

Ym = + Rm) ^; yRci = (1 + -Rhoi + ^hci-Rhcii) ^; yHcii = -Rhoi(1 + -Rhci + -Rhoi-Rhcii) ^, (16) 

where 

R.^^, (17) 

Ti is the photoionization rate of species i G {H i. He i. He 11}, 

du^a.iiy), (18) 

and ai is the (case A) recombination coefficient to all atomic levels of species i. The recombination rate of the next 
ionization state i + 1 (e.g., if i is H i then i -I- 1 is H 11) is neUi^iai, where the electron number density rig is 

rie = 71h(1 - Yai) + ?^Ho>HcII + 2nHc(l - Y^i^j - IhcIi)- (19) 

In the case of a highly ionized medium with i?Hi, ^^HcI, -RhoII 3> 1, the densities of He i and He 11 can be expressed in 
terms of the H i density as 

'^Hci ^ riHc Rm j.2Q^ 
nni ?^H ^HcI^HcII 

and 

"HcII ^ "-He Rm ^2]^-) 

nm n-H i?Hcii 

For optically thin systems, the above relations with = J^, clearly give the ratio between the column densities of 
different ions. Note how the quantity 

?y = iVHcii/iVHi (22) 
is independent on gas density only as long as the optically thin approximation holds, while the ratio 

C = Nnei/Nm (23) 

is always density dependent. 
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4.1. Slab approximation and fitting formulae 

An iterative solution to the equations of cosmological radiative transfer that included a detailed numerical calculation 
of the ionization and temperature structure of individual absorbers at every timestep would be a very computing- 
intensive task. To properly treat the self-shielding of LyC radiation, in Paper II we modeled absorbers as semi-infinite 
slabs, developed a "steplike" approximation to the function /^(iVHi), and used an analytical escape probability formalism 
to include the recombination emission from absorbers. Fardal, Giroux, & ShuU (1998) solved the local radiative transfer 
problem via an integral equation (the Milne solution for a gray atmosphere) for the number of photoionizations at any 
optical depth in a given slab. They also devised an approximation formula to rj that closely followed the numerical 
results. Faucher-Giguere et al. (2009) have recently generalized the treatment of Fardal et al. (1998) and applied a 
similar fitting formula to the results of a code that self-consistently solves the photoionization equilibrium balance, 
including the influence of recombination radiation. Here, we follow a similar method: under the assumption of Jeans 
length thickness for the absorbers, we solve the ionization and thermal structure of a slab of finite width illuminated 
by an external isotropic radiation field J^, and derive analytical approximations for the ratios rj and C as a function 
of iVfji. Details of our calculations are provided in the Appendix. We parameterize the external background flux as 
a power- law, J„ = Jgi2{i' / 1^912)^°' , and as in Faucher-Giguere et al. (2009) divide the intensity above 54.4 eV by a 
factor of 10 to mimic a cosmological UV filtered spectrum. 

Figure H shows the resulting ratios 7?(iVHi) and C(-^hi) for a range of input spectra and for the representative intensity 
value at 1 ryd of 10~^^ ergscm"^ s~^ Hz"^ sr~^. The function 77 remains constant at low H i columns, as long as the 
optically thin approximation holds. As the neutral hydrogen column increases, the slab first becomes optically thick 
to He ii-ionizing radiation, and rj increases. Slabs with even larger columns become optically thick to H i LyC: they 
are characterized by a highly ionized surface layer and an almost fully neutral core. This is the reason for the rapid 
decrease of r] after the peak, and the consequent trend of ( toward the neutral limit, C "-Hc/'^h- As in Fardal et al. 
(1998) and Faucher-Giguere et al. (2009), we calculate the column iVnoii from the equation 

riRc 7-912 p _ ^ , ^228 p . X 

-; 1 I A ^HI - T228 + -, , p -KHell, (,^4j 

4nH 1 + Atqi2 1 + Bt28S 

where T912 = A'^hiO'912, ^228 = A^HcIiO'228, and A and B are constants fitted to our numerical results. To make use of the 
above expression, one must further specify the ionization rates F^ to be used in the Ri terms together with a relation 
between electron density and A^hi- It is in this second step that our approach differs from that of Faucher-Giguere et 
al. (2009). These authors used the optically thin limit for r^. which provides a poor approximation to the numerical 
results. Here, we first compute the ionization rates in the optically thin limit, 

4-TT I 

dv^aM. (25) 
hv 



^thin 



and derive rf^™ using equation (pl|). For a given (input value) of thi, we then calculate t^'J.j". We then write a 
first-order approximation to the He 11 ionization rate at the face of the slab, 

r?fen =\dv l^e-Hon(^) ^^^„(,). (26) 



I di.^e-^-MaHi(i^). (27) 



The analogous expression for H i is 



Finally, we compute the Ri factors for H i and He 11 in equation (g4|) as 

R. = 0-5rf'" + 0-5rf- (23^ 

neat 

The recombination rates depend on the gas temperature. We found that our numerical results can be fit with the 
simple scaling 

T = (2 X 10^ K) J9V2,-22, (29) 

where J9i2,-22 = ^912/10^^^ ergs cm~^ s~^ Hz"-'^ sr~^, is adequate for our purposes. The weak dependence on J is 
related to the fact that, for log A^hi ^ 16, cooling is largely provided by coUisionally-excited line radiation rather than 
by recombinations. With this simplified treatment, we have been able to fit the numerically obtained values of rj for 
a broad range of input spectra. The best-fit curves shown in the left panel of Figure |^ have been obtained taking 
A = 0.02 and B = 0.25 in equation (|2|), and 

ne = 3.0 X 10- W3 (A^Hi,i7.2)'/' (r^jt-12)'^' (30) 
for the electron density. Here, AHia7.2 = A^m/lO^'^-^ cm'^ and r^^in = r^hin/j^Q-12 g-i^ r^j^g ^^^^^ relation can 
be derived assuming Jeans length thickness for the absorbers and optically thin photoionization equilibrium (Schaye 
2001; see also Faucher-Giguere et al. 2009). A simple approximation for ( can be also derived, once 7j is obtained. We 
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use equation (pO|), ( = rj/Rj^ci, with FhoI = for iVni ^ 10^^ cm ^. At larger columns we apply a linear (in log 

space) extrapolation to the limiting vale C — !> riHo/^-H assumed to be reached at A'hi = 10^^ cm~^. 
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Fig. 3. — Left panel: The ratio r] = A''Hcn/-'VHi as a function of A^hi at rcdshift 3. The illuminating spectrum has intensity = 
■^912 (i'/i'9i2)~" , with J9i2,-22 = 1 and, from bottom to top, spectral slopes o = 0, 1, 2, with a break of a factor of 10 at 54-A-eV. Solid 
curves: full numerical photoionization calculations. Dashed curves: our analytical approximations for r] based on equation (124). Dotted 
curves: optically thin limit. Right panel: same for f = A''hci/^hi- 



5. RECOMBINATION EMISSIVITY 

In § 2 we have seen how background photons absorbed through a Lyman series resonance cause a radiative cascade 
that ultimately terminates either in a Lya photon or in two-photon 2s ^ Is continuum decay. In this section we use 
the detailed photoionization structure of absorbing systems to calculate the reprocessing of background LyC radiation 
by the clumpy IGM via atomic recombination processes. We include recombinations from the continuum to the ground 
state of H i, He i, and He ii, as well as He ii Balmer, two-photon, and Lya emission. Using the formalism developed 
in the Appendix, the recombination flux at the slab surface. 



(31) 



can be written as 



= t; dxju{x) I dfic 
^ Jo Jo 



-T„{x)lp. 



The emission coefficient from a generic recombination process. 



dx2v(x)E2{Ty{x)). 



(32) 



2v{x) = hv(j>^arne{x)ni+i{x) = hiy(j>^ — ni{x)Ti{x), 

a. 



(33) 



where (/)^ is the normalized emission profile and ar is the relevant recombination coefficient, is proportional to the 
density of species i, times the rate at which it absorbs ionizing photons (L^), times the fraction of recombinations that 
lead to the radiative transition under consideration (the ratio ar/ai). The emission profile of free-bound recombination 
radiation can be computed via the Milne detailed-balance relation, which relates the velocity-dependent recombination 
cross section to the photoionization cross section, while a delta-function profile is sufficient for bound-bound transitions. 
The cosmological proper recombination emissivity for the relevant recombination process can then be computed by 
integrating over the distribution of absorbers, 



e,iz) = 2\dz/cdt\ / dNm /(^hi, ^)^^.(^hi), 



(34) 
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where the factor 2 accounts for the two surfaces of a slab. Using equations (|3^) and (p3|), and denoting with Ni = 
riidx the species i column density of the absorber, the recombination emissivity becomes 

^ foo fNi(NHi) 

e,{z) = \dz/cdt\h,^(l),— dNHif{NHi,z) dN[T,{N'^)E2{T,{N'^)). (35) 

"i Jo Jo 

As with the ionization and thermal structure of individual absorbers, it is not practical to perform a self-consistent, 
iterative, numerical evaluation of the recombination emissivity at every timestep in the cosmological code. To derive 
a simple analytical formula to the emergent radiation from absorbers, we make use of the fact the number of ionizing 
incident photons that are absorbed saturates in the optically thick regime, and approximate the second integral on 
the rhs of equation (|3|) as (cf. Fauchcr-Giguere et al. 2009) 

I{Ni) = J dN[ T{Nl)E2{T^{Nl)) « (O.Srf + 0.5rf Nt (l - e"^'/^^) . (36) 

Here Nt is the column density of ion i above which the recombination emission saturates. As shown in Figure 3, 
the above formula works especially well in the case of LyC recombination re-emission from H i and He ii, where self- 
absorption by the emitting ion dominates the local reprocessing of recombination radiation. Our best-fit parameters 
to the full numerical results for H i. He i, and He ii LyC recombinations are Nt — 6.5 x 10^^ (i'/i'9i2)^ '^ cm~^, 
Nt = 1.2 X 10^^ ( I'/ 1^504)"'^'^ cm~^, and Nt = 2.3 x 10"'^^ (j^/j^228)^'^ cm~^, respectively. (In the case of non-ionizing H i 
recombination Lya and two-photon emission, we find Nt = 6.5 x 10^^ cm~^.) 

The emergent recombination flux from He ii BalC, two-photon, and Lya depends on the helium (emission) as well as 
hydrogen (absorption) ionization structure. With the adopted column density distribution, however, recombinations 
into He ii are dominated by absorbers in the range of columns 10^^ < A^hi ^ 10^^ cm in these systems, H i 
absorption can be neglected and a simple approximation can be found by setting Nt = 2.3 x 10^^ cm^^. The fit 



at large H i columns is actually improved by multiplying the rhs of equation (36) by exp[— min(rgi2, 1.3)](i'gi2/j^)*''^. 
A comparison between the results of the full numerical integration of the localradiative transfer equation and our 
analytical approximations to the recombination radiation from He ii BalC, two-photon, and Lya are shown in the 
right panel of Figure ^. Note that, in our calculations, we have again assumed that He ii Lya photons diffuse into the 
wings and then escape subject only to continuum absorption. 




Fig. 4. — Left panel: LyC recombination radiation from quasar absorbers as a function of Afjji- The points depict the.results of the full 
numerical integration of the local radiative transfer equation, while the lines show our analytical approximations (eq. pq ). Green squares: 
H I LyC at 912 A. Blue triangleAi-Rc II LyC at 228 A. Red circles: He I LyC at 504 A. The quantity plotted is the integral I{Ni) (cq. pd) 
multiplied by Niiif{Niii, z) (cq. |l2[ ) at ^ = 3, showing the contribution of optically thin and optically thick absorbers to the LyC emissivify. 
For this comparison we assumeoan illuminating spectrum with Jqi2,-22 = 1 and spectral slope a = 1, with a break of a factor of 10 at 
54.4 eV. Right panel: same as left panel, but for recombination re-emission from He II BalC at 13.6 eV [blue squares), Lyo at 40.8 eV 
(green triangles), and two- photon continuum at 20.4 eV (red circles). 
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6. QUASAR UV EMISSIVITY 

The only sources of ionizing radiation included in CUBA are star-forming galaxies and quasars. For the quasar 
comoving emissivity at 1 ryd, egi2{z)/{l + z)^, we use the function 



£912(^) 

(1 + ^)3 



(10^^-^ ergss-^Mpc-3Hz"i)(l + z) 



4.68 



exp(-0.28z) 
cxp(1.77z) + 26.3' 



(37) 



which closely fits the results of Hopkins, Richards, & Hcrnquist (2007) in the redshift interval 1 < z < 5. The UV 
SED is given by the broken power-law 

^ ri.-o-44 (A > 1300 A); 

^-«j^-i.57 (A < 1300 A) ^-^^^ 
(Vanden Berk et al. 2001; Teller et al. 2002). 
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Fig. 5. — Quasar comoving emissivity at 2 keV (dashed line) and 10 kcV [solid line) in units of 10^"^ crgss ""^ Mpc ''Hz ^. latter 
has been computed following Ueda et al. (2003) and Silverman et al. (2008), the former using the procedure outlined in Section [7.2[ 



7. QUASAR X-RAY EMISSIVITY 



The extrapolation of the steep UV power-law in equation (^8|) to higher energies is unable to reproduce the X-ray 
properties of the quasar population as a whole, as recorded in the cosmic X-ray background (XRB). The XRB may play 
a unique role in regulating the thermodynamics and ionization degree of intergalactic absorbers. In a photoionized 
IGM, soft X-rays between 0.5 and 0.9 keV are responsible for the highest ionization states of metals like carbon, 
nitrogen, and oxygen. At early epochs. X-rays penetrate regions that are optically thick to UV radiation, providing a 
source of heating and ionization. They could make the IGM warm and weakly ionized prior to the era of reionization 
breakthrough (e.g.. Oh 2001; Venkatesan, Giroux, & ShuU 2001; Madau et al. 2004; Ricotti & Ostriker 2004; Kuhlen 
& Madau 2005). Compton scattering of hard XRB photons may be a source of heating for highly ionized low-density 
intergalactic gas (Madau & Efstathiou 1999). 

Deep X-ray surveys aided by optical identification programs have shown that the bulk of the XRB is produced 
by a mixture of unobscured "Type 1" and obscured "Type 2" AGNs (Mushotzky et al. 2000; Giacconi et al. 2001), 
as predicted by XRB synthesis models constructed within the framework of AGN unification schemes (e.g., Setti & 
Woltjer 1989; Madau, Ghisellini, & Fabian 1994; Comastri et al. 1995; Gilh, Comastri, & Hasinger 2007). Here, we 
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compute the total X-ray emissivity from Type 1 and Type 2 AGNs following a modern version of the original approach 
by Madau et al. (1994). 

7.1. Intrinsic hard X-ray luminosity function 

According to Ueda et al. (2003), who combined various surveys from the HEAO 1, ASCA and Chandra satellites, 
the hard 2-10 keV quasar luminosity function (HXLF) follows a luminosity-dependent density evolution with a cutoff 
redshift (above which the evolution stops) that increases with luminosity. At the present epoch, the intrinsic (i.e., 
before absorption) HXLF of all AGNs (including both Type I's and Type 2's) is best represented by 



in the luminosity range lO''^'^ — 10*^'^ ergss ^, where = 2190 Gpc ^ and = lO'*'^-^^ ergss ^. This changes with 
cosmic time (for redshift up to 3) as 

</)(L,z) = (/.(L,0)e(z,L), (40) 



where the evolution factor is 



and 




x^i-h.A^^ (41) 



(L > La) 
r(i/La)"-33'5 {L<La) 



= (i) = r.*../r ^0.335 )rZ ri' (42) 



Here, ei = 4.23, 62 = -1.5, z* = 1.9, and La = lO'*'^-^ ergss'^ (Ueda et al. 2003). An extension of the HXLF up to 
z 5 by Silverman et al. (2008) shows a steeper decline in the number of z > 3 AGNs with an evolution rate similar 
to that found by studies of optically-selected QSOs. The new fit requires a much stronger evolution above the cutoff 
redshift, 62 = —3.27, than previously found by Ueda et al. (2003, 62 = —1.5). In the following, we shall use Ueda et 
al. (2003) HXLF best fit parameters together with the Silverman et al. (2008) value for 62. 
For the intrinsic spectrum before absorption, we assume the standard power-law multiplied by an exponential, 

5b cx-B-^expf-^V (43) 



with a = 0.9 (Nandra & Pounds 1994). The high-energy cutoff, E,. = 460 keV, is fixed by the shape of the XRB 
turnover above 30 kcV. These seed photons are then reflected towards the observer by a semi-infinite cold disk close 
to the primary emitter. This reflection component, commonly detected in the X-ray spectra of nearby Seyfert galaxies 
(Nandra & Pounds 1994), is comparable to the direct flux around 30 keV, decreases rapidly towards lower energies, 
and flattens the overall spectral slope above 10 keV (Lightman & White 1988). 

7.2. AGN emissivity after absorption 

According to the AGN unification scheme, obscuring matter at a distance of several parsecs from the central pow- 
erhouse blocks our line of sight to the active nucleus. When our view is unobscured, we see a Type 1 AGN; when our 
view is occulted, photons of all energies from the far IR to several keV are absorbed, and in these bands we can only 
detect the nucleus in scattered light. Ueda et al. (2003) found the following expression for the observed (normalized) 
distribution of absorbing iVfj columns: 

r2 - (5 + 2e)/(l + e)ip (20.0 < logiVn < 20.5), 
/(i,iVH) = h/(l + e)V (20.5 < logiVn < 23.0), (44) 

[e/{l + e)ip (23.0 < logiVH < 24), 

where the parameter 

V'(i) = min{^max, max[V'44 - ^(logL - 44), 0]} (45) 
accounts for the fact that the fraction of absorbed sources is smaller at higher luminosities. Here, V-'max = (l+£)/(3 + e), 
e = 1.7, ■044 = 0.47, /3 = 0.1, and J^^ f{L, N^) dlogNu = 1. Sources absorbed by a column larger (smaller) than 
10^^ cm^^ are defined as X-ray Type 2 (Type 1) AGNs. It is assumed that "Compton-thick" AGNs with columns 
logA^H > 24 are not present in samples detected below 10 keV. Such a population is added by extrapolating the A'h 
function above logA^^f > 24, keeping the same normalization up to logA^u — 25 as well as the same cosmological 
evolution of Compton-thin AGNs (Ueda et al. 2003). 

We then follow Madau, Ghiscllini, & Fabian (1993) and model the thick blocking material that covers most of the solid 
angle around the central X-ray source as a homogeneous spherical cloud of cold material and column Nn. The radiation 
transfer is computed with a Monte Carlo code constructed using the photon-escape weighing method of Pozdnyakov, 
Sobol', & Sunyaev (1983). We set the electron temperature equal to zero, use the full Klein-Nishina scattering 
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cross section, adopt the bound-free opacity associated with standard cosmic-abundance material from Morrison & 
McCammon (1983), and ignore the iron Ka emission hne in the spectra. Each Monte Carlo run uses 10^ photons, 
and produces as output an absorbed spectrum, SEiN-a). After being reprocessed by cold material along the line of 
sight, the emergent specific intensity forms a hump, whose position and width arc determined by the competition of 
bound-free absorption at low energies, and Compton downscattering and exponential roll-off of the primary spectrum 
at high energies (Madau et al. 1993). A small spectral component, equal to 2.5% of the primary incident power 
and representing the flux scattered into the line of sight by electrons in the warm ionized medium, is added to the 
transmitted Type 2 flux. The absorbed spectra are then averaged over the A'n-distribution corresponding to a given 
luminosity, and normalized to the unabsorbed 2-10 keV flux. 



1^0 S E{Nii)f{L,Nii)d\ og iVn 

l2-10keV ^EdE 



Se{L) = "V"^^^ ' ^ (46) 



to yield a flux normalized, luminosity-dependent, average AGN SED. The X-ray proper emissivity as a function of 
redshift is then obtained by simply integrating over the HXLF, 

eE{z) ^ (1 + zf "^'^^ SE{LmL,z)LdL, (47) 

where we set Lmin = lO'^^-^ ergss"^, and imax = 10*® ergss"^. 

The model described above is able to reproduce a number of X-ray observations, from the evolution of AGNs in the 
soft and hard X-ray bands, to the XRB. The quasar comoving emissivity at 2 keV and 10 kcV is plotted in Figure |^, 
while a global fit to the XRB is shown in the left panel of Figure ^. The absolute XRB fiux is still affected by rather 
large uncertainties: our model reproduces well the background intensity measured by HEAO-1 and BeppoSAX, but 
the HEAO-1 A2 data arc lower by about 20% with respect to the determinations by, e.g., XMM and RXTE at energies 
below 10 keV. Figure ^ (right panel) depicts the broadband quasar comoving emissivity per logarithmic bandwidth, 
-\- zY , as a function of photon energy from the optical to hard X-rays. In terms of energy output, the composite 
spectrum for A < 5000 A is characterized by two broad bumps, one in the UV at 10 eV and another in the X-ray 
region at 30 keV (a third peak in the mid- infrared, see, e.g., Sazonov, Ostriker, & Sunyaev 2004, can be neglected 
for the present purposes). While X-rays dominate the energy ouput at z = 0, the peak of the emitted power moves 
increasingly towards the UV at redshifts above 1. 




Fig. 6. — Left panel: The cosmic XRB spectrum and the predicted contribution from AGNs. Grey points: HEAO-1 A2 HED data (Gruber 
et ah 1999). Dark green points: HEAO-1 A4 LED (Gruber et al. 1999). Cyan points: Rossi-XTE (Revnivtsev et al. 2003). Blue point: 
0.25 keV soft XRB intensity from iJOSAT shadowing experiments (Warwick & Roberts 1998). Red bowtie: HEAO-1 A4 MED (Kinzer et 
al. 1997). Blue bowtie: RO S AT PSPC data (Georgantopoulos et al. 1996). Light green bowtie: BeppoSAX (Vecchi et al. 1999). Purple 
and yellow bowties: Newton-XMM (Lumb et al. 2002; E)e Luca & Molendi 2004). Solid line: our synthesis model spectrum, produced 
by a mixture of absorbed (logA'^jj > 22, short-dashed line) and unabsorbed (logA^H < 22, long-dashed line) AGNs. See text for details. 
Right panel: The broadband quasar comoving emissivity per logarithmic bandwidth, vey/{l -\- z)^ (in units of 10"^^ ergs Mpc~^), as a 
function of photon energy E from the optical to hard X-rays. The composite spectrum is shown at redshifts 2 = 0, 1, 3, 5. 
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8. GALAXY EMISSIVITY 

Star-forming galaxies arc expected to play a dominant role as sources of hydrogen-ionizing radiation at z > 3 as 
the quasar population declines with lookback time. To compute the LyC emissivity from galaxies at all epochs, we 
start with an empirical determination of the star formation history of the universe following Madau et al. (1996). We 
adopt the far-UV (FUV, 1500 A) luminosity functions of Schiminovich et al. (2005) in the redshift range < z < 2, 
of Reddy & Steidel (2009) at z = 2.3 and 3.05, and of Bouwens et al. (2011) at redshifts 3.8, 5.0, 5.9, 6.8, and 8.0. 
All were integrated down to L,„in = 0.01 using Scliechter function fits with parameters (0*,L*,a) to compute the 
dust-reddened galaxy FUV luminosity density pfuv0, 

/•CXD 

PFuy{z)= L(l){L,z)dL = T{2 + a,0.0l)<p^L^. (48) 

Jo. OIL, 

Here a denotes the faint-end slope of the Schechter parameterization and T is the incomplete gamma function. We 
used a = -1.6 at < z < 2, a = -1.73 at z = 2.3 and z = 3.05, and a = -1.73,-1.66,-1.74,-2.01,-1.91 at 
z = 3.8,5.0,5.9,6.8,8.0, respectively (sec Schiminovich et al. 2005; Reddy & Stcidcl 2009; Bouwens et al. 2011). 
Dust attenuation was treated using a Calzetti et al. (2000) extinction law, with the function 

measuring the magnitudes of attenuation suffered at frequency ly and redshift z. For the luminosity- weighted obscu- 
ration at 1500 A we take 

... [1 (0 < z < 2); 

AFUv(^)-|2.5iog[(i + i.5/(z-l)] (z>2). ^^^> 

The above expression reproduces at 2 < 2 the FUV "minimum dust correction factor" of 2.5 from Schiminovich et al. 
(2005), the dust correction factors of 2.38±0.59 and 2.0±0.62 at z = 2.3 and z = 3.05 from Reddy & Steidel (2009), and 
the decreasing dust attenuation at higher redshift from Bouwens et al. (2011). The dust-corrected luminosity densities 
were smoothed with an approximating function and then compared with the results of spectral population synthesis 
models as follows. The GALAXEV library of Bruzual & Chariot (2003) provides the age-luminosity evolution for a 
simple stellar population (SSP) at different wavelengths. The FUV luminosity density (before dust obscuration) at 
time < of a "cosmic stellar population" characterized by a star formation rate density SFRD(t) and a mctal-cnrichmcnt 
law Z(t) is given by the convolution integral 

PFUv(i)= / SFRD{t -T)lFvv[T,Z{t-T)]dT, (51) 
Jo 

where Zfuv['''j Z{t — r)] is specific luminosity radiated at 1500 A per imit initial stellar mass by an SSP at age r and 
metallicity Z{t — t). We use SSPs of decreasing metallicitics with redshift according to 

Z{z) = Z0lO-°-i5" (52) 

(Kcwley & Kobulnicky 2007), for a Salpcter IMF between 0.1 and 100 Mq. Starting from an initial guess, the function 
SFRD(t) was adjusted in an iterative fashion until the computed FUV luminosity densities as a function of redshift 
provided a good match to the data. The best-fitting star formation history, 

SFRD(.) = ''^"^^l^J™" Mpc-^ (53) 

is shown in Figure |^ (left panel), together with the observed luminosity densities adopted in this study. The latter 
have been converted to ongoing star formation rate densities according to 

SFRD(t) = /C X pFuv(0; ^ = 1-05 x IQ-^^, (54) 

where pFUV is expressed in units of ergss"^ Mpc"'^ Hz"-'^ and SFRD is in units of Mq yr~^ Mpc~'^. This approximate 
transformation makes use of the basic property that the FUV continuum in galaxies is dominated by short-lived 
massive stars, and is therefore a direct measure, for a given IMF and dust content, of the instantaneous star formation 
rate. The conversion factor K. in the equation above reproduces to within 2% the results of the synthesis models above 
redshift 2 given the adopted star formation and metal enrichment history (cf. Madau, Pozzetti, & Dickinson 1998). 
At redshift < z < 1, decreases from 1.16 x 10~^* to 1.10 x 10~^^. Note that these newly derived conversion factors 
are between 21% and 33% smaller than the widely used value, /C = 1.4 x 10~^*, quoted by Kennicutt (1998) (and 
based on the calibration by Madau et al. 1998), the differences rcfiecting updated stellar population synthesis models 
and subsolar metallicitics at high redshifts. 

^ In this section we use tlic notation Pv{z) = e,y{z)/{l + z)^, i.e., the term luminosity density is synonymous with comoving specific 
emissivity. 
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Fig. 7. — Left: The cosmic history of star formation. The FUV data points from Schiminovich et al. (2005) [blue dots), Reddy 
& Steidel (2009) {red squares) and Bouwens et al. (2011) {magenta pentagons) have been converted to instantaneous star formation 
rate density using the conversion factor K. = 1.05 X 10~^* (see text for details). The best-fitting star formation history, SFRD(2) = 
[6.9 X 10~^ + 0. 14(2/2. 2)l'5]/[l + (2/2.7)*-^] M0 yr-l Mpc'^, is plotted with the solid blue curve. Right: Comoving galaxy emissivity (in 
units of 10^'^ ergs s"-' Mpc"'' Hz"-'^) of 1 ryd photons cscapixig into the IGM {dashed line), for an escape fraction (/esc) = 1.8xlO~*(l+2)'^ ''. 
The solid line shows the best-fit QSO emissivity of eq. (pTI) for comparison, while the dot-dashed line shows the total quasars -|- galaxies 
emissivity. 



Once the star formation history has been determined, we use stehar synthesis models to compute the dust-reddened 
frequency-dependent UV emissivity as 

p^{t)=:C{t) I ^YY{D{t-T%[T,Z{t-T)]dT. (55) 
Jo 

We take C{t) = lO~°-^"^(^'*) at all photon energies below 1 ryd, and C{i) = (/esc) above the Lyman limit. In our 
treatment, the escape fraction (/esc) is a free parameter that incorporates local continuum absorption by hydrogen, 
helium, and dust. It is the angle-averaged, absorption cross section- weighted, and luminosity- weighted fraction of 
ionizing radiation that leaks into the IGM from star-forming galaxies: the escaping radiation is produced not by 
sources in a semiopaque medium but by a small fraction of essentially unobscured sources (e.g., Gnedin, Kravtsov, & 
Ghen 2008). In the "minimal reionization model" discussed in detail in the next section, the escape fraction of photons 
between 1 and 4 ryd is assumed to be a steeply rising function of redshift (see also Inoue, Iwata, & Deharveng 2006), 

(/esc)- 1.8 X 10-4(1 + z)3-4, (56) 

and is zero above 4 ryd. The expression above yields an escape fraction at z ~ 3.3 of 2.6%, comparable to the recent 
upper limit for L > Lyman break galaxies of Boutsia et al. (2011). The relatively low values of (/esc) implied 
by the above expression in the redshift interval from z = 2 (0.8%) to 2; = 5 (8%) are dictated in our model by the 
need to reproduce the hydrogen- ionization rates inferred from flux decrement measurements (see Fig. |^ below). In 
the same redshift range, the escape fraction of ionizing radiation from star-forming galaxies hosting a 7-ray burst is 
measured to be (/esc) < 7.5% (95% c.l.) (Ghen, Prochaska, & Gnedin 2007), in agreement with our expression. The 
high values predicted by equation (M) above redshift 7, in excess of 20%, are needed to compensate for the decline 
in the star formation rate density Etna to reionize the IGM at early enough epochs. The resulting galaxy emissivity 
of 1 ryd photons escaping into the IGM is shown in the right panel of Figure 0. Galaxies dominate over QSOs at 
all redshifts z > 4, and make a negligible contribution to the ionizing background at z < 3. The total comoving 
emissivity from quasars -I- galaxies decreases only weakly from z = 3 to z = 5, and is fairly flat afterwards. This 
trend is consistent with the conclusions reached by Bolton & Haehnclt (2007) and Fauchcr-Giguere et al. (2008a) from 
empirical measurements of the Lya forest opacity. 

8.1. Lya emission from galaxies 

Stellar population synthesis codes do not typically include nebular line emission. Here, we provide a simple estimate 
of the Lya emission from hydrogen recombinations in the interstellar medium of galaxies. In case B recombination, 
about 68% of all the absorbed LyG photons will be converted locally into Lya (Osterbrock 1989). The Lya proper 
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volume emissivity can then be written as 

ea{z) = hiya6{iy - Va)na{z), (57) 

where 

/•CO T 

n<,(z) = 0.68(1 -(/e.c)) / j^e^iz) (58) 

and ei, is the proper volume emissivity from galaxies. We assume here that Lya suffers the same dust extinction 
as LyC, a simple treatment that is unlikely to capture the complex radiative transfer physics of the Lya line as it 
propagates through the dusty ISM (see, e.g., Caplan & Deharveng 1986; Neufeld 1991; Dijkstra 2009; Scarlata et al. 
2009; Dayal, Ferrara & Saro 2010). Inserting equation ( |57| ) into (|^) yields the additional flux observed at Vo < t'a from 
galaxy Lya : 

where 1 + Za ~ {i'a/i^o){^ + Zq)- We have neglected collisionally excited Lya emission, as this is only about 10-20% of 
the recombination term (Dayal et al. 2010). A similar contribution is also expected in the emitted spectrum of dense 
absorbers like the SLLSs and DLAs, while coUisional excitation is always negligible in lower column density systems. 




1 2 3 4 5 6 7 0123456789 
redshift redshift 

Fig. 8. — Left: The hydrogen photoionization rate, Fhi, from ^ = 1 to 2: = 7. Solid curve: quasars + galaxies model. The dashed 
curves depict the individual contributions of the QSO population {blue) that dominates at low redshift and of the galaxy population 
(red) that reionize the IGM at early times. Circles: empirical measurements from the Lya forest effective opacity by Bolton & Hachnelt 
(2007). Triangles: same by Becker, Ranch, &; Sargent (2007) (their lognormal model). Squares: same by Faucher-Giguero et al. (2008a). 
Pentagons: same using the quasar proximity effect by Calverley et al. (2011). Right: The hydrogen photoheating rate per ion, 'Hhi {upper 
set of curves), and the He n photoheating rate, WhcH {lower set of curves), from the present epoch to z = 9. AH photoheating rates are 
expressed in units of 10~^^ eV s~^. Solid lines: quasars + galaxies. Dashed lines: quasar-only. The addition of a galaxy component boosts 
the H I rate and decreases the He n rate. The dotted line shows the Compton heating rate per electron in units of 10~^* eV s^^. 



9. BASIC RESULTS 

This section gives a quick overview of the main results generated by the upgraded CUBA radiative transfer code, 
using the formalism and parameters described above. CUBA solves the radiative transfer equation (^ by iteration, as 
its right-hand term implicitly contains J in the recombination emissivity and in the effective helium opacity. 

9.1. Photoionization and photoheating rates 

The total optically thin photoionization rate of hydrogen, Fhi, is shown in Figure || as a function of redshift (left 
panel). For comparison, we have also plotted the individual contributions of the QSO population that dominates at low 
redshift and of the galaxy population that reionize the IGM at early times, together with the empirical measurements 
from the Lya forest effective opacity by Bolton & Haehnelt (2007), Becker, Rauch, & Sargent (2007), and Faucher- 
Giguere et al. (2008a), and from the quasar proximity effect by Calverley et al. (2011). The fractional recombination 
contribution to Fhi increases from 9% at 2; = to 18% at z = 4 to up to 37% at z > 7: it does so because the mean 
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free path of recombination photons decreases with lookback time and a smaller fraction of such photons gets redshifted 
below the ionization threshold before capture (Faucher-Giguere et al. 2009). While the total H i photoionization rate 
provides a good match to the data, we note that there are large systematic uncertainties in the measurements as these 
depend on the assumed IGM temperature and gas density distribution. 

The right panel of the same figure depicts the optically thin photoheating rates per ion (see eq. A7 with J^^, = for 
the definition) of hydrogen, Hhi, and He ii, HhcII, for the quasars + galaxies and quasar-only models. The addition 
of a galaxy component boosts the H i rate as it increases the emissivity of hydrogen-ionizing photons at fixed H i 
opacity (the latter being determined by the observations). The opposite is true for He ii photoheating (as well as 
He II photoionization), as galaxies do not contribute to the emissivity above 4 ryd: this increases the predicted He ii 
opacity (again at fixed H i opacity) and causes a large break in the background spectrum at 4 ryd and a smaller 
photoheating rate. While the Compton heating rate per electron is many orders of magnitude (about 7 dex at redshift 
3) below the H i photoheating rate (note the different normalization of the heating rates plotted in Fig. it is a 
non-negligible source of heating for very underdense, highly ionized regions: the Compton heating rate for intergalactic 
gas at ovcrdensity 0.1, temperature T = 10'' K, and redshifts z = (1, 2, 3) is (53, 16, 4)% of the total photoheating rate. 
Table 3 tabulates the optically thin photoionization and photoheating rates of hydrogen and helium predicted by our 
"quasars -|- galaxies" model for use, e.g., in cosmological hydrodynamics simulations of the Lya forest. 

9.2. Background spectral energy distribution 

Figure H shows the spectrum of the radiation background as a function of redshift for a "quasar-only" model, together 
with the old results from Paper H. The new spectra are characterized by a lower UV flux (by as much as a factor of 
3 at 1 ryd and z = 3), smaller spectral breaks at 1 and 4 ryd because of the reduced H i and He ii LyC absorption, a 
sawtooth modulation by the Lyman series of H i and He ii that becomes more and more pronounced with increasing 
redshift, and a flatter soft X-ray spectrum. 

The addition of radiation from galaxies has little effect on the ionizing background at redshifts below 3, as shown 
in Figure |l^, a consequence of our adopted redshift-dependcnt escape fraction. At higher redshifts the impact is more 
dramatic: a large boost at 1 ryd is associated with a much sharper He ii sawtooth and He ii absorption edge. As 
noted above regarding the photoheating rates, this arises because galaxy spectra arc truncated at 4 ryd, and the large 
increase in the H-ionizing emissivity from the early galaxy population is not accompanied by a similar increase at the 
He II edge. The net effect is a larger He ii opacity at fixed H i opacity. At z > 5, the He i opacity of the IGM also starts 
building up (it is negligible at lower redshifts), and a small He i absorption edge can be discerned in the spectrum of 
the background at 24.6 eV. At z > 3, the sawtooth modulation produced by resonant absorption in the Lyman series 
of intergalactic He ii (see Fig. |ll|) is clearly a sensitive probe of the nature of the sources that keep the IGM ionized, 
and may be a crucial ingredients in the modelling of the abundances of metal absorption systems (Madau & Haardt 
2009). The analogous sawtooth modulation produced by the H i Lyman series becomes significant above redshift 6 
(see Fig. Oh, and may affect the photodissociation of molecular hydrogen during cosmological reionization (Haiman, 
Rees, & Loeb 1997). 

Figure |l3| compares the broadband spectrum of the total extragalactic background light (EBL) from quasars and 
galaxies, predicted by CUBA at z = 0, with current EBL observations from the mid-IR to the 7-rays. 

9.3. A "minimal reionization model" 

It is interesting at this stage to use the quasar and galaxy ionizing emissivities of § 6, 7, and 8 and track the evolution 
of the volume filling factors of ionized hydrogen and doubly ionized helium regions in the universe as a function of 
cosmic time. As shown in Paper HI, the volume filling factor of H 11 regions, Qhii, is equal at any given instant t 
to the integral over cosmic time of the number ionizing photons emitted per hydrogen atom by all radiation sources 
present at earlier epochs, 

minus the number of radiative recombinations per ionized hydrogen atom, 
Here 

?^io„(^)= / (./csc)T^e.(t) (62) 

with (/osc) = 1 hi the case of quasars, (nn) = 1-9 x 10^^(1 + z)^ cm^'^ is the mean hydrogen density of the expanding 
IGM, and (tree) is the volume-averaged hydrogen recombination timescale, 

{tree) = [x{nH)aBC]-\ (63) 

where as is the recombination coefficient to the excited states of hydrogen, x = 1-08 accounts for the presence of 
photoelectrons from singly ionized helium, and Cigm = ('^hii)/(^hii)^ is the clumping factor of ionized hydrogen. 
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The cosmic background photoionization and photoheating rates. 
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Fig. 9. — The broadband spectrum of a "quasar-only" cosmic background between 5 A and 5,000 A at epochs z = 0,1, 3, and 5. 
The new models {black curves) are compared with the old results of Paper II {blue curves). The intensity is expressed in units of 



IQ- 



ergscm s Hz sr . The vertical thin linos indicate the positions of the H I and He II Lya and Lyman limit. 



Differentiation yields the H i "rcionization equation" of Paper III, 

C^QhII _ "-ion <3hII 

dt (nn) {tree) ' 

and its equivalent for expanding He iii regions, 

^QhcIII ?^ion,4 QhoIII 



dt 



(riHe) (^rec,Ho} ' 



(64) 



(65) 



where /iion,4 now includes only photons above 4 ryd (which are mostly absorbed by He ii), and the recombination 
timescale of doubly ionized helium, (ircc.Hc): is the about 6 times shorter than the hydrogen recombination timescale 
if H II and He iii have similar clumping factors. We will not attempt here to model the rcionization of He i, as this 
occurs nearly simultaneously to and cannot be readily decoupled from that of H i. The rcionization equation equation: 
1) describes the transition from a neutral universe to a fully ionized one in a statistical way, independently, for a given 
emissivity, of the emission histories of individual radiation sources; 2) assumes that the mean free path of ionizing 
photons is much smaller than the horizon, i.e., that they are absorbed before being redshifted below the ionization 
edge; and 3) includes in the source term only those photons above the Lyman limit that escape into the IGM ((/esc) — 1 
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Fig. 10. — The broadband spectrum of a "quasars + galaxies" cosmic background at redsliifts z = 1.1,3,4.9, and 5.9 [black 
curves). The new "quasar-only" model of Fig. m is plotted for comparison {red curves). The intensity .Jv is expressed in units of 
10~^^ crgscm"^ Hz~^ sr~^. The vertical thin lines indicate the positions of the H I and He n Lya and Lyman limit. 

in the case of quasars) . Photons that are absorbed in loco by dense interstellar gas do not enter in the source term, nor 
does the interstellar absorbing material contribute to the recombination rate. The volume-weighted clumping factor 
reflects only the nonuniformity of the ionized low-density IGM, the repository of most of the baryons in the universe, 
and its use in the recombination timescale is justified when the size of the ionized regions is large compared to the 
scale of the clumping. 

When Qhii 1 (the "pre-overlap" stage), individual ionization fronts propagate from star- forming early galaxies 
into the low-density IGM. The neutral phase shrinks as Qhii grows and H ii regions start to overlap. The radiation 
field remains highly inhomogeneous until the reionization process is completed at the "overlap epoch", Qhii = 1, when 
all the low-density IGM becomes highly ionized. Pockets of neutral gas remain in collapsed systems during the entire 
"post-overlap" stage (Gncdin 2000) and may manifest themselves as the SLLSs or DLA systems in quasar absorption 
spectra. We have integrated equation (^) assuming a gas temperature of 2 x 10** K and a clumping factor for the 
intergalactic medium of 

CiGM = l + 43z-i-^\ (66) 

This is equal to the expression for Cioo (the clumping factor of gas below a threshold ovcrdensity of 100) found at 
z > 6 in a suite of cosmological hydrodynamical simulations by Pawlik, Schayc, and van Schcrpenzeel (2009). These 
authors found that photoionization heating by a uniform UV background greatly reduces clumping as it smoothes out 
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Fig. 11. — A zoom-in of the "quasars + galaxies" (black curves) and "quasar-only" (red curves) cosmic background spectrum at redshifts 
z = 3, 3.4, 4.1, and 4.6 showing the sawtooth modulation of the mctagalactic flux between 220 and 320 A produced by resonant absorption 
in the Lyman series of intcrgalactic He II. The vertical thin lines indicate the positions of the He II hya and Lyman limit. 



small-scale density fluctuations, and that the clumping factor at z = 6 is insensitive to the redshift at which the UV 
background is actually turned on (as long as reheating occurs at > 9). We use an overdensity of 100 to differentiate 
between dense gas belonging to virialized halos and the diffuse intcrgalactic gas, and assume that the collapsed mass 
fraction is small. We also extrapolate equation (|6^) down to 2 > 2, and assume the same clumping factor for H 11 and 
He III. 

The results of this "minimal reionization model" are shown in Figure |lj. Cosmological H 11 regions driven by 
star- forming galaxies overlap at redshift 6.7, and the hydrogen in the universe is half- ionized (by volume) at redshift 
10. He III regions driven by quasars overlap much later, at redshift 2.8, and their filling factor is only 4% at redshift 
5. These overlap epochs are consistent with the SDSS spectra of 2; ~ 6 quasars (Fan et al. 2006b), with numerical 
simulations of H i reionization (Gnedin & Fan 2006), and with observations of the He 11 Lya forest at z < 3 (see, e.g., 
Worseck ct al. 2011; ShuU ct al. 2010; Fechner et al. 2006; Heap et al. 2000; and references therein). A simple probe 
of the reionization history is the integrated optical depth to electron scattering Tes, which depends on the path length 
through ionized gas along the line of sight to the CMB as 

r dz' 

Tesiz)^aTC . [QhIi(»h) + QueuinHe) + 2QHeIII (^He)] (67) 

JQ J^[i- + Z ) 
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Fig. 12. — A zoom-in of the "quasars + galaxies" cosmic background spectrum at redshifts z = 5.2, 6.1, 7.2, and 8.0 showing the sawtooth 
modulation of the metagalactic flux between and 890 and 1300 A produced by resonant absorption in the Lyman series of intergalactic H I. 
The vertical thin lines indicate the positions of the H I Ly« and Lyman Hmit. The dashed line shows the same spectrum without sawtooth 
for comparison. 



(Wyithe & Locb 2003), where ut is the Thomson cross section. The seven-year PFMylPresults imply Tos = 0.088±0.015 
(Jarosik et al. 2011). Our minimal reionization model assumes QhoII = Qhii and yields an electron scattering opacity 
to the epoch of reionization of Tcs ~ 0.084, in good agreement with the observations. 

The outcome of our minimal reionization model is rather sensitive to the assumed escape fraction of hydrogen- 
ionizing radiation at early epochs; this exceeds 50% at z > 9 and reaches unity at z = 11.6. Had we assumed a 
maximum /esc of 50% instead, the same model would yield Qhi = 1 at z = 6.2 and Tcs = 0.06. We also remark that, 
in the pre-overlap era, the background spectra shown in Figures |[ |o[ ^ and |2| have only a formal meaning, as 
they describe a space-averaged radiation field that is in reality highly inhomogeneous. Recent spectra taken by the 
Cosmic Origins Spectrograph on the Hubble Space Telescope exhibit patchy He ii Gunn-Peterson absorption, with a 
mean He ii/H i abundance ratio that is 47 ± 42 at 2.4 < z < 2.73, and 209 ± 281 at z > 2.73 (ShuU et al. 2010). In 
the rcdshift interval 2.4 < z < 2.73, our background spectrum yields a He ii/H i abundance ratio (in the optically thin 
limit) around 50-70, in good agreement with the observations. The predicted mean He ii/H i ratio increases rapidly 
towards high redshift, to (280,494,887,1615) at z = (3.42,3.64,3.87,4.1) as galaxies start dominating the ionizing 
emissivity and the spectrum of the UVB steepens. The evolution of the He ii abundance and the fluctuating spectrum 
of the cosmic UVB in the pre-overlap era will be the subject of a subsequent paper. 
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-4-3-2-10 1 2 3 
log hi/ (keV) 

Fig. 13. — The predicted broadband cxtragalactic background light, vj^, from quasars and galaxies at 2 = 0, compared with empirical 
determinations at different wavelengths. Red points: the optical-near IR EBL from HST and ground-based galaxy counts (Madau & 
Pozzctti 2000). Blue points: tho-|mid-IR EBL from IRAC-Spitxer galaxy counts (Fazio et al. 2004). The X-ray data points are explained 
in details in the caption of Fig. B. 



9.4. The UVB: uncertainties 

The background spectra computed in the previous section are sensitive to a number of poorly determined input 
parameters. In this section we brieflv discuss just a few of the uncertainties inherent in our synthesis modelhng of 
the UVB. The left panel of Figure |l5| shows the adopted comoving quasar emissivity at 1 Ryd (eq. together 
with the determinations by Meiksin (2005), Cowie, Barger, & Trouille (2009), Bongiorno et al. (2007), Willott et al. 
(2010), and Siana et al. (2008). The poorly known faint-end slope of the quasar luminosity function at high redshift, 
incompleteness corrections, as well as the uncertain spectral energy distribution (SED) in the UV, all contribute to 
the large apparent discrepancies between different measurements. 

The escape fraction of ionizing radiation that leaks into the IGM and its unknown redshift evolution is the major 
source of uncertainty in the determination of the galaxy contribution to the UVB. To better gauge the impact of 
this parameter on our synthesis model, we show in the right panel of Figure ^ the spectrum of the comoving galaxy 
emissivity at four different epochs. The relatively large leakage of LyC photons assumed at early times is a crucial 
ingredient of our "minimal reionization model" , which yields an optical depth to Thomson scattering in agreement with 
WMAP results. A smaller escape fraction at high redshifts would lead to too-late reionization, while a significantly 
larger escape fraction at lower redshifts would produce a hydrogen photoionization rate that appears to be too high 
compared to the observations (see Fig. ||). 

We have also checked that, for a given IMF, uncertainties in the stellar population synthesis technique are relatively 
small. Figure ^ shows the emission rate of hydrogen-ionizing photons for an SSP, calculated as a funtion of age with 
the GALAXEV models of Bruzual & Chariot (2003), the Starburst99 models of Leitherer et al. (1999), and the FSPS 
models of Conroy, Gunn, & White (2009). While the three packages use different stellar evolution tracks and spectral 
libraries, the total number of ionizing photons emitted agrees to within 10%. The figure also illustrates the significant 
effect of stellar metallicity: an SSP of metallicity 1/50 of solar emits 60% more hydrogen- ionizing photons over its 
Hfetime than a solar metallicity SSP (Salpeter IMF, GALAXEV package). 

Finally, we address the effect of a change in the effective opacity of the IGM. In our parameterization, the shape 
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redshift redshift 



Fig. 14. — Evolution of the H II (left panel) and He III {right panel) filling factors as a function of redshft for our "minimal reionization 
model" (see text for details). The dotted curve in the left panel depicts the cumulative electron scattering optical depth of the universe, in 
units of 10%, as a function of redshift. 




redshift E (eV) 



Fig. 15. — Left: Comoving specific emissivity at 1 ryd (in units of 10^"^ ergss"^ Mpc"'' Hz~^) measured from different quasar surveys. 
Solid line: the best-fit function in eq. (^). Dashed line: Cowie et al. (2009). Dot-dashed line: Bongiorno et al. (2007), using a 912 A to 
4400 A flux ratio of 0.31. Filled circles: Meiksin (2005), PLE modeh Empy circles: Meiksin (2005), PDE model. Empty star: Siana et ah 
(2008). Filled square: Willott et al. (2010) at z = 6. Right: Comoving galaxy emissivity per logarithmic bandwridth (in units of 10'^'' ergs 
s~^ Mpc~'^) escaping into the IGM, as a function of photon energy, at four different redshifts. Note the large, time-evolving break at 13.6 
eV. 



of the /(-/Vhijz) distribution over the column density range of the LLSs is adjusted at every redshift for continuity 
with the SLLSs. As detailed in § 3, this procedure yields the slopes (3 — 0.47, 0.61, 0.72, 0.82 at redshifts 2 = 2,3, 4, 5, 
respectively. To gauge how an uncertainty in /(A^hi) translates into an uncertainty in the UVB, we have run CUBA 
with the fixed values of /? = 0.2 and /? = 1 in the column density interval 10^^'^ < A^hi < 10^^ cm~^. The resulting 
UVB at z = 3 is shown in the right panel of Figure ^ Compared to our fiducial model, the fiat /3 = 0.2 distribution 
generates a hydrogen photoionization rate (Thi) that is 29% lower, and a He ii photoionization rate (FhoIi) that is 
22% higher. This is because a larger opacity at 1 Ryd from the LSSs results in a harder background spectrum, which 
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in turn produces a smaller He ii opacity at 4 Ryd. Conversely, in the steep /3 = 1 case, Fhi increases by 8% and FhoII 
decreases by 10%. Notice, however, that the former model would significantly underestimate the 1 Ryd photon mean 
free path compared to the measurements of Prochaska et al. (2009). 




Fig. 16. — Left: Emission rate of hydrogen-ionizing photons for a SSP of total mass 1 Mq, Salpeter IMF, and solar metallicity, as a 
function of age. Solid line: results from the GALAXEV package of Bruzual & Chariot (2003). Long-dashed line: same for Starburst99 
(Leitherer et al. 1999). Short-dashed line: same for the FSPS package of Conroy et al. (2009). Dotted line: same as the solid line, but for a 
metallicity 1/50 of solar. Right: Uncertainties in the broadband spectrum of the "quasars + galaxies" UVB at redshift 3. Black line: our 
fiducial model. Green line: a model with a slope of /3 = 0.2 in the column density distribution of LLSs. Red line: same for /3 = 1.0. The 
intensity is expressed in units of lO^^'^ ergsem^'^ s^^ Hz^^ sr^^. 



10. SUMMARY 

In this paper we have presented improved synthesis models of the evolving spectrum of the UV/X-ray diffuse back- 
ground, updating and extending our previous results. Five new main components have been added to our cosmological 
radiative transfer code CUBA and discusses in details: (1) the sawtooth modulation of the background intensity from 
resonant line absorption in the Lyman series of cosmic hydrogen and helium; (2) the X-ray emission from the obscured 
and unobscured quasars that gives origin to the X-ray background; (3) a piecewise parameterization of the distribution 
in redshift and column density of intergalactic absorbers that fits recent measurements of the mean free path of 1 ryd 
photons; (4) an accurate treatment of the absorber photoionization structure, which enters in the calculation of the 
helium continuum opacity and recombination emissivity; and (5) the UV emission from star-forming galaxies at all 
redshifts. The full implications of our new population synthesis models for the thermodynamics and ionization state 
of the Lya forest and metal absorbers will be addressed in a subsequent paper. Here we have provided tables of the 
predicted H i and He ii photoionization and photoheating rates for use, e.g., in cosmological hydrodynamics simu- 
lations of the Lya forest, a new mctallicity-dependent calibration to the UV luminosity density-star formation rate 
density relation, and presented a "minimal cosmic reionization model" in which the galaxy UV emissivity traces recent 
determinations of the cosmic history of star formation, the luminosity-weighted escape fraction of hydrogen-ionizing 
radiation increases rapidly with lookback time, the clumping factor of the high-redshift intergalactic medium follows 
recent determinations of hydrodynamic simulations that include the effect of photoionization heating, and Population 
HI stars and miniquasars make a negligible contribution to the metagalactic flux. The model has been shown to provide 
a good fit to the hydrogen-ionization rates inferred from flux decrement and quasar proximity effect measurements, to 
predict that cosmological H ii (He iii) regions overlap at redshift 6.7 (2.8), and to yield an optical depth to Thomson 
scattering. Tea = 0.084 that is agreement with WMAP results. 

Our new background intensities and spectra are sensitive to a number of poorly determined input parameters and 
suffer from various degeneracies. Their predictive power should be constantly tested against new observations. We 
are therefore making our redshift-depende nt UV/X emissivities and CUBA outputs freely available for public use at 
jittp : //www ■ ucolick ■ org/~pmadau/ CUBA . 
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APPENDIX 

THE IONIZATION AND THERMAL STATE OF INTERGALACTIC ABSORBERS 

In this Appendix we describe the numerical calculations of the ionization structure of individual absorbers outlined 
in § ^ We approximate each absorber as a semi-infinite slab with uniform hydrogen density tih and thickness equal 
to its Jeans length, 

L = ^'K^kTl{Gpinn^) = ^mV^cn^'^ tI'^ (Al) 

(Schaye 2001; Fauchcr-Gigucre ct al. 2009), where the adiabatic index and mean molecular weight for a monoatomic, 
fully ionized gas of primordial composition are 7 = 5/3 and \i = 0.59, respectively, and Ti^ is the gas temperature in 
units of 10* K. The gas mass fraction is set at its universal value, jg = fib/f^M = 0.16. Each slab model provides, 
given the gas density tih and an external isotropic radiation field J^, the columns iVui, A^HcI, and iVucii. The slab 




Fig. 17. — The slab approximation for individual absorbers: a plane-parallel slab of thickness L illuminated by an isotropic external 
radiation field J. 



geometry reduces the local radiative transfer to a 2-D problem (see Fig. n^. Setting /i = cos 6*, where B is the angle 
between the incident ray and the normal at the slab, the specific intensity7j/(a;, [i) at coordinate x within the slab is 
given by the following implicit solution to the radiative transfer equation: 

r/,(0,Ai)e--'^(-)/^ +3i^/J=dx'>(x')e-[-"(-)--''(-')l/M (0 < < ^2); 
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New synthesis models of the cosmic UV/X-ray background 



Here, j^, is the recombination radiation emission coefficient along the ray, 

T^{x) ^y2a-i{v) / dx'ni{x'), (A3) 

is the optical depth normal to the slab, and the sum is to be taken over all relevant species. Integration over solid 
angle yields the local specific brightness Ju{x) 

M^) = ^ J dm, =:|{iJ2[r,(x)]+ii;2[r,(L)-T,(.T)]} + i^^rfx'^^^i?i(|T,(.T)-r,(xO|), (A4) 

where Ei and E2 are the first and second order exponential integral functions, respectively. The symmetry of the 
problem allows us to study just half of the slab, as for all ions n(x) = n(L — x). With this in mind, we can rewrite the 
external contribution to J,{x) (the first term on the right hand side of eq. |A4D as 

^ {E2[tAx)] + i?2[2r,(L/2) - r,(.T)]}, (A5) 



(0 < a; < i/2), and the contribution of recombinations (the second term on the right hand side of eq. A4) as 

1 /'^^'dx':^{i?i(|r,(x)-r,(x')|) +i?i(2r,(i/2)-T,(:E')-T.(x))}. (A6) 

2 Jo 47r 

Thermal equilibrium in the absorbers is obtained by balancing photoionization heating with a number of different 
energy loss mechanisms: free-free, coUisional excitations and ionizations, recombinations, Compton cooling against the 
cosmic microwave backgroud, and adiabatic cooling from cosmic expansion. The photoheating rate (per ion) is given 
by 

fdiy ^h{iy - i^M^), (A7) 



hi/ 

where hui is the ionization potential of species i. The equation of thermal equilibrium can then be written as 

nrnHm + ^Hei'HHei + mdi'HHcU = ?^eAc + ^kTHln-a + 4nHc)/Ai, (AS) 

where A,, is the atomic cooling rate and the last term accounts for adiabatic cooling. Our treatment of the thermal state 
of cosmological absorbers is over-simplified in two aspects: 1) the lowest density IGM is not in thermal equilibrium; 
and 2) dense quasar absorbers do not expand with the Hubble flow. We find, however, that: 1) adiabatic cooling 
is never important for absorbers with A^hi ^ 10^^ cm~^, so the inclusion of this process does not alter the thermal 
balance in this column density regime; and 2) in low the density, optically thin IGM, the functions 77 = A'hcIi/A'hi 
and (to a lesser extent) C, = Nuei/Nm are independent on gas temperature, and hence the assumption of thermal 
equilbrium should not affect our results. 

The ionization and thermal structure of the slab, for a given J, and nn, are solved by iteration. We start by assuming 
an almost fully ionized slab at T = 10^ K, and compute the ionization rates Ti setting the recombination emissivity 
to zero. We then solve for the ionization structure within the slab, compute updated ion fractions and opacities, and 
evaluate new ionization rates including recombination radiation. The process is iterated until convergence to better 
than 0.1% is obtained at every point of the slab. At each iteration the size of the slab is changed accordingly to the 
Jeans criterion (eq. Al). 



